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Abstract 

We develop a set of equations to describe the population dynamics of many interacting species 
in food webs. Predator-prey interactions are non-linear, and are based on ratio-dependent 
functional responses. The equations account for competition for resources between members 
of the same species, and between members of different species. Predators divide their to- 
tal hunting/foraging effort between the available prey species according to an evolutionarily 
stable strategy (ESS). The ESS foraging behaviour does not correspond to the predictions 
of optimal foraging theory. We use the population dynamics equations in simulations of the 
Webworld model of evolving ecosystems. New species are added to an existing food web 
due to speciation events, whilst species become extinct due to coevolution and competition. 
We study the dynamics of species-diversity in Webworld on a macro-evolutionary timescale. 
Coevolutionary interactions are strong enough to cause continuous overturn of species, in 
contrast to our previous Webworld simulations with simpler population dynamics. Although 
there are significant fluctuations in species diversity because of speciation and extinction. 
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very large scale extinction avalanches appear to be absent from the dynamics, and we find no 
evidence for self-organised criticality. 

1 Introduction 

Understanding coevolution within communities of interacting species is one of the greatest 
challenges in the study of ecological systems. Two different sets of issues are involved in mod- 
elling such communities. Questions regarding food web structure, the nature of predator-prey 
interactions, competition for resources, and population dynamics apply on an ecological time 
scale comparable to the lifetime of individual organisms. Questions regarding evolutionary 
change of species, introduction of new species to the food web by speciation processes, and 
removal of species due to extinction apply on an evolutionary time scale orders of magnitude 
longer than the lifetime of an organism. We argue that these two sets of questions are nev- 
ertheless related and that they need to be considered within the same framework. In order 
to understand food web structures we need to understand the way in which the diversity of 
organisms in the ecosystem evolved. In order to understand the way the set of species in a 
food web will coevolve we need to understand the nature of the competitive interactions and 
predator-prey relationships between them. 

The Webworld model, introduced by Caldarelli, Higgs & McKane (1998), and studied 
further here, is an attempt to model the two timescales simultaneously. The model considers 
a set of species, each of which has a set of morphological and behavioural features that 
determine the way it interacts with all the other species, and hence the positions of links in 
the food web. Population dynamics equations are used to determine the way the population 
sizes of all the species change over ecological time. In one evolutionary time step of the model 
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one new species is added to the food web, and the populations of all the species then change 
in response. The new species sometimes adds stably to the ecosystem, sometimes dies out due 
to competition with existing species, and sometimes causes the extinction of other species. 
The diversity of species within the ecosystem thus changes on the evolutionary time scale due 
to speciation and extinction. 

In our previous paper (Caldarelli et al (1998), henceforward Paper I), we considered the 
properties of the food webs generated by Webworld, including the number of basal, interme- 
diate and top species in the web, the number of links per species, and the number of trophic 
levels. These properties have been measured in real food webs (e.g. Cohen et al., 1990; Hall 
& RaffaeUi, 1991; Goldwasser & Roughgarden, 1993; Martinez & Lawton, 1995), and thus an 
extensive amount of ecological data exists with which we were able to compare the results of 
the model. The Webworld model generates food webs with properties that are in reasonable 
agreement with those of real webs, given the large uncertainties inherent in measurements on 
real webs. The model also makes predictions about the way food web properties will change 
as a function of ecological parameters such as the rate of input of external resources to the 
ecosystem, the efficiency of transfer of resources from prey to consumer at each level of the 
food chain, and the strength of competition between species for the same resources. As such, 
we feel that the model goes considerably further than other theoretical models of food web 
structure, such as the cascade model (Cohen, 1990; Cohen et al, 1990), which are based on 
constructing random graphs. 

The Webworld model also makes predictions about the dynamics of species diversity that 
can be compared with the evidence from the fossil record. There has been considerable interest 
in macro-evolutionary models recently, generated by the claim that extinction dynamics are 
related to the concept of self-organised criticahty (Bak & Sneppen, 1993; Sole et al. 1996; 
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Sole et al, 1997; Amaral & Meyer, 1999). The idea is that the avalanches of extinctions 
visible in the fossil record can be expected to arise from the internal coevolutionary dynamics 
of the system, and thus one does not need to postulate catastrophic external events such 
as meteorite strikes or climate changes in order to explain the extinctions. In simulations 
of the Webworld model in Paper I we found that large extinction avalanches could occur in 
situations when the ecosystem was poorly adapted to the external conditions, but that as 
time went by extinction events became smaller and rarer. The ecosystem tended towards a 
'frozen' food web of mutually well-adapted species that could not be invaded by new species. 
These results therefore did not support the idea of criticality. 

Whilst models such as that of Bak & Sneppen (1993) have the merit of being (deliberately) 
very simple, we have found that the dynamics of ecosystem models depends substantially on 
the way that such models are set up, and we feel that it is important to attempt to include some 
degree of realism in the models if one wishes to draw conclusions about the real world. One 
of the aims of this Paper is to develop a general set of equations for population dynamics that 
deals with competition between species and predator-prcy interactions in a food web which 
can have any arbitrary structure of links generated by the evolutionary process. The equations 
used here are based on ratio-dependent functional responses (Arditi & Ginsburg, 1989; Arditi 
& Akgakaya, 1990; Arditi & Michalski, 1995) and represent a considerable improvement on 
those used in Paper I in the way they treat competition between predators for a given prey. 
Another important change is that increased adaptation of predators leads to a decrease in 
prey population size in the current model, whereas this was not so in the model used in Paper 
I. This leads to a continuous overturn of new species replacing old ones, in contrast to the 
frozen state found in Paper 1. Although the stationary state is now a dynamical one, we still 
find no evidence for self-organised criticality. 
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The outline of the paper is as follows. In section 2 we define the Webworld model as 
studied both here and in Paper I. In section 3 we present the new equations for population 
dynamics. We show that these equations satisfy several logical requirements of general food 
web models. We also give an interpretation of the choice of diet of predators in terms on 
evolutionary stable strategies. Section 4 gives details of the simulations of the Webworld 
model using the new equations. In sections 5 and 6 we investigate the long term evolutionary 
dynamics, firstly in the absence of predators and then in webs with predation. We conclude 
in section 7 with a discussion of the implications of these changes and with a comparison with 
other evolutionary models. 



2 The Webworld model 

This section describes the basis of the Webworld model introduced in I. Details that differ from 
Paper 1 will be mentioned specifically. Our model is a stochastic one, since the characteristics 
of the species and the speciation events are chosen randomly, however we use deterministic 
dynamical equations for the population sizes of each species. A species is defined by the set of 
its morphological or behavioral characteristics. We construct a species in the model by picking 
L features out of a pool of K possible features. These features represent morphological and 
behavioural characteristics, which might be, for example, "nocturnal", "having sharp teeth" 
and "ability to run fast" , however, we do not assign particular biological attributes to each 
feature: they are just integers which run from 1 to K. In our simulations we take L = 10 and 
K = 500 for illustrative purposes. 

The matrix of scores mapi describes how useful one feature, a, is against any other feature, 
j3. The K X K matrix rriap is antisymmetric (i.e., rUais = —mpa) and is taken to consist 
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of random Gaussian variables with mean zero and unit variance. These are chosen at the 
beginning of a simulation run and do not change during that particular run. The score Sij of 
one species i against another species j is then defined as 



where a runs over all the features of species i and (3 runs over all the features of species j. 
Thus the score of one species against another is essentially just the sum of the scores of the 
relevant features against each other. A positive score Sij indicates that species i is adapted for 
predation on species j, whilst a zero score means that there is no predator- prey relationship 
between the species. The scores will be used in the equations for population dynamics de- 
scribed in the next section. The external environment is represented as an additional species 
which is assigned a set of L features randomly at the beginning of a run and which does not 
change. Species having a positive score against the external environment represent primary 
producers of the ecosystem. 

The model consists of a changing set of species that may feed on each other and on external 
resources. External resources are input at a constant rate R and are distributed amongst 
species as a function of their scores, in a way that is discussed later. These resources are then 
tied up in the ecosystem as potential "food" in the form of prey for predator species. For 
simplicity, we measure resources and population sizes in the same units, so that Ni{t) denotes 
the number of individuals for species i at time t or alternatively the amount of resources 
invested in species i at this time. 

The short time dynamics is described by a set of equations giving the change of population 
size of any one species in terms of the population sizes of the other species in the ecosystem. 
The form of these equations is to be discussed in the next section. A new species is created 




(1) 
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by a speciation event from one of the existing species. This is carried out by choosing a 
parent species at random, and introducing a new daughter species into the ecosystem that 
differs from the parent species by one randomly chosen feature. The new species begins with 
a population size of 1, and 1 is subtracted from the population of the parent species. The 
populations of all the species are determined by iterating the population dynamics equations. 
If the population size of any species falls below one, that species is removed from the system, 
and so rendered extinct. The population dynamics simulation is continued until all surviving 
species reach an equilibrium population size, or until a defined large time period is reached 
without the populations having reached equilibrium. This completes one evolutionary time 
step of the model, and the program proceeds to add another new species. In order to prevent 
multiple copies of identical species from arising, each time a new species is added, a check is 
carried out to ensure that the set of features of the new species is not already represented in 
the ecosystem. 

A minor change between the simulations here and those in Paper I is the way that species 
are chosen to undergo speciation: in the original model they were chosen with a probability 
proportional to their population size, here they are chosen randomly. As we will discuss below, 
this does not lead to qualitative changes in the behaviour of the model. The major change in 
the model is the form chosen for the population dynamics. This is discussed in detail in the 
next section. 

3 Population dynamics 

We wish to develop a set of population dynamics equations which is general enough to deal 
with any food web structure. There have been many models of population dynamics that 
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discuss only two or three species — e.g. plant plus herbivore plus carnivore, or two consumer 
species competing for the same resource. However many of these models are not easy to 
generalise to the multiple species case. Most species in an ecosystem are both predators and 
prey and are in competition with several other species. We require equations which include 
all these effects at the same time. 

Let the rate at which one individual of species i consumes individuals of species j be 
denoted by Qijit). This is usually called the 'functional response', and it depends in general 
on the population sizes. We suppose that the population size of each species satisfies an 
equation of the form: 



The first term on the right represents a constant rate of death of individuals in absence 
of interaction with other species. The final term is the sum of the rates of predation on 
species i by all other species, and the middle term is the rate of increase of species i due to 
predation on other species. Where there is no predator-prey relationship between the species 
the corresponding rate gij is zero. For primary producers the middle term includes a non-zero 
rate of feeding on the external resources. The factor A is less than 1, and is known as the 
ecological efficiency. It represents the fraction of the resources of the prey that are converted 
into resources of the predator at each stage of the food chain. Throughout this paper, we 
have taken A = 0.1, a value accepted by many ecologists (Pimm, 1982). We have deliberately 
chosen the form of Eq. (^) to be the same for all species. We do not want to define different 
equations for primary producers, herbivores, and carnivores etc, because species can change 
their position in the ecosystem as it evolves, and most species are both predators and prey. 
For simplicity, we have set the death rate to be equal for all species and the value of A to be 



dN,{t) 
dt 



N,{t) + \Y^N,g,,{t)-Y.N,g,,{t). 
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equal for all species. In a more complex model we could have allowed these quantities to be 
functions of the sets of features of each species and then these parameters would have been 
subject to evolution in the same way that the interactions scores between species are subject 
to evolution. The choice of the death rate to be unity in (§) essentially sets the time scale 
for the population dynamics: the time appearing in this equation has been scaled so that the 
coefficient of —Ni{t) is one. 

Equation (0) is different from Eq. (5) in Paper I, which had been designed to be as simple 
as possible, and to have only one stationary state. First, Eq. (5) in Paper I is discrete in time, 
while Eq. is continuous. However, this difference is of minor importance, as Eq. (|^) has to 
be discretized anyway for performing computer simulations, giving Eq. (|16|) below, which has 
a similar form to Eq. (5) of paper I in the case of large time steps. At = 1. The second and 
main difference to Paper I consists in the fact that the term describing decrease in population 
size due to predation in Paper I was chosen to be independent of the rate at which the species 
is consumed by its predators. Also, the form of the functional response in the consumption 
term had not been chosen according to ecological considerations. 

The main question to be addressed in this section is how to choose a reasonable function 
Qij that is applicable for all the situations that arise in the ecosystem. Since the final form 
we use is relatively complex, we will describe several particular cases first and build up to the 
general case. 

For a single predator i feeding on a single prey j we suppose 

This is known as a ratio-dependent functional response (Arditi & Ginsburg, 1989; Arditi & 
Michalsi, 1995), because Qij can be written as a function of the ratio of prey to predators 



if both top and bottom are divided by the predator population Ni. When the prey is very 
abundant, gij = Sij/b, i.e. each predator feeds at a constant maximum rate. When predators 
are numerous compared to the available prey, there is competition between predators, and the 
rate at which each individual predator can feed on the prey becomes limited by the amount 
of prey. In this limit the combined rate of consumption of all predators is NiQij = Nj. This 
situation is known as donor control. Arditi & Akgakaya (1990) have shown that interference 
between predators is significant and that the ratio-dependent functional response can be 
applied to a wide range of real species. 

In our model, the same equation is used to treat consumption of external resources by 
a primary producer. In this case i is the primary producer, and the external resources are 
j = 0, with a value of A^o that is kept fixed and equal to R/ A. By writing down the differential 
equation (2) for the case of a single primary producer species 1, we find that the population 
size A^i reaches a stationary value when Xqiq = 1. Hence, the equilibrium population size is 

N, = (A5io - fe)A^o/^io , (4) 

provided this is positive, i.e. species 1 can only survive if Sio > b/\. Thus b is an important 
parameter of the model that determines the minimum score necessary for a consumer to 
survive on a given resource. With the choice of A^o given above, A^i tends to R if it is very 
well adapted {Sio ^ 1). The parameter R represents the fixed rate of supply of non-biological 
resources, principally sunlight. These resources are renewed constantly and hence are never 
depleted. Also they cannot accumulate if not used, hence there is no differential equation for 
the rate of change of A"o. 

If there are additional species competing with i for predation on species j, equation (§) 
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can be generalised as follows: 

^''^^^ ^ bNj{t)+Ekt!sk,N,{ty 
The sum in the denominator is over all species k which prey on j, including i itself, i.e. it 
is over all species for which Skj > 0. The additional predator populations are present in the 
denominator because each individual of species i is in competition with the other species as 
well as with other members of its own species. The factor is introduced to represent the fact 
that competition between members of the same species for a resource is usually stronger than 
competition between members of different species. Thus aki < 1 when k and i are not equal 
and ttfcfc = 1 for all k. We will in addition suppose that — (^ik- Although addition of this 
extra factor complicates the equations, it is actually essential in order to permit coexistence 
of competing species. As an example consider two species 1 and 2 competing for external 
resources. In this case: 

' bNo{t) + SM{t) + auS2oN2{t) ' ^ ^ 

/,N 'S'2oA^o(t) 

^''^ > ~ bNoit) + SMt) + auS,oN,{t) ' ^ ' 

In the stationary state Xgio = 1 and Xg2o = 1, hence 

No{X{Sw - ai2-S2o) - b{l - a^)) 
5io(l - 0:12) 



A^o(A(5'2o - Qi2'S'io) - 6 (1 - ai2)) 
520(1 - ai2) 



ivov^V^20 — <'«i2'Jio; — "v-^ ~ "12;; /n\ 

^TTi — 7r\ ■ y^) 



For coexistence of the two species both the above must be positive, thus the species can only 
coexist if 

- (1 - a,2){Sw - b/X) < S20 - ^10 < (1 - ai2)(^2o - b/X). (10) 
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Therefore the range of the difference between scores for which coexistence is possible is pro- 
portional to 1 — ai2- If the competition between species is reduced {au is reduced) it becomes 
easier for species to coexist on the same resources. If ai2 = 1 only species with identical scores 
can coexist. Since in general there is more than one predator per prey in real food webs, it 
is necessary to introduce the a parameter into the model. The result that species can only 
coexist if between species competition is weaker than within species competition is also found 
in other models (e.g. Renshaw (1991) pp 137-139). 

For the purpose of our simulations we will suppose that the strength of competition de- 
pends only on the degree of similarity between the species: 

aij = c+ {1 - c)qij, (11) 

where c is a constant such that < c < 1, and with qij being the fraction of features of species 
i that are also possessed by species j. 

We now consider the case of a single predator with more than one prey. It might seem 
that we could use equation (|^) for each prey j. However this is unsatisfactory, as explained 
below. In fact we use: 

bN,it) + s,,umity ^ ' 

where we have introduced the factor fij, which is the fraction of its effort (or available searching 
time) that species i puts into preying on species j. These efforts must satisfy J2j fij = 1 for 
all i. The importance of introducing the efforts can be understood by considering a single 
predator i = 3 with two prey j = 1 and j = 2 of population sizes A^^i and In the 

particular case where the prey are equivalent from the predator's point of view (i.e. 6*31 = 532), 
the dynamics of the predator population should be identical to the case where there is just 
one prey population of size A^^i + A^2- This is a 'common sense' condition that has been 
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emphasised by Arditi & Michalski (1995) and Berryman et al (1995), who have shown that 
many dynamical equations used previously did not satisfy the condition. In our case, since 
the prey are equivalent, the predator sets its efforts to be proportional to the population sizes: 
/31 = Ni/{Ni + N2) and /32 = N2/{Ni + N2). Calculating the predation rates from equation 
(12) the total input to the predator can be shown to be 



which is the same as for a combined species of population A^^i + A^2- If the efforts had not been 
introduced, i.e. equation (|) had been used instead of equation ([T^), this condition would not 
be satisfied. 

There is an additional reason why it is necessary to introduce the fij. The rate of decrease 
of a prey population j caused by a predator i is NiQij. If the simple equation (|^) is used, then 
the effect of the predators on the prey is very large when the predator population is large 
compared to the prey. When a new species evolves, it always begins from small population 
numbers. Usually, there is an existing species which has some ability to act as a predator 
on the new species. The existing predator has a relatively large population because it must 
already be successfully feeding on an established prey species in the ecosystem. Thus the new 
prey species suffers from an enormous level of predation and almost always becomes extinct 
as soon as it is introduced, even if it is substantially better adapted than the established prey. 
This prevents a diverse ecosystem from evolving. This problem is solved by introducing the 
efforts, since initially the predator puts very little effort into feeding on the new prey because 
there are so few of them. The effect of the predator on the new prey is thus in proportion 
to the prey's population size. This permits newly-evolved species to enter the ecosystem in a 
reasonable way. 
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We now require a rule by which predators assign their efforts to different prey when the 
prey are not equivalent. We suppose that the efforts of any species i are chosen so that the 
gain per unit effort Qijl fij is equal for all prey j. If this were not true, the predator could 
increase its energy intake by putting more effort into a prey with higher gain per unit effort. 
This choice of efforts leads to the condition 



Ek9ik{t)' 



m = (14) 



It is shown in the appendix that this choice is an evolutionarily stable strategy (ESS) (Parker 
& Maynard Smith, 1990). If the population has efforts chosen in this way, there is no other 
choice of efforts that can do better, i.e. no other strategy than can become more common 
if it is rare. When prey are equivalent, the ESS solution reduces to setting the efforts in 
proportion to the prey population sizes, as above. 

Combining all the above considerations we arrive at the following general form for the 
functional response that is used in the Webworld simulations in this paper: 

a (f) = s^^hmAt) .... 
bN,it) + EkC^HSk,fkjm,ity ^ ' 



with the efforts given by Eq. ([1^). 

We have previously shown that a restricted form of Eq. (|T5|) was invariant under aggre- 
gation of identical prey. It is straightforward to demonstrate that this holds for the general 
form (^) too. The invariance under aggregation of identical predators can also be shown. If 
predator i and predator I are identical, we have Sij = Sij and aij = aij for all j, and an = 1. 
The combined effect of the two species on prey j is therefore 



NiQij + NiQij 



bNj + Sij{Nifij + Nifij) + J2kyti,i o^kifkjSkjNk 
14 



which is obviously identical to the effect of one predator species of population size Ni + A^"; 
and effort 

{NJi^ + Nifi,)/{N, + Ni). 

Therefore our equations Eq. and Eq. ([T5|) satisfy the logical requirements of invariance 
under aggregation of identical species. We now go on to consider the practicalities of imple- 
menting these equations in simulations. 



4 Implementation of the Webworld program 

Solution of time-dependent differential equations involves a numerical algorithm such as the 
Runge-Kutta method which integrates forward by small time steps. We require to simulate the 
dynamics of large numbers of species over large times, hence the efficiency of the algorithm is 
important. To speed up the computer simulations, we used a discrete version of the dynamics. 



Ni{t + At) = N,{t){l - At) + At 



\J2N,g,,it)-J2N,g,,it) 



(16) 



with a time-step At = 0.2, which is quite large. The discrete version of the dynamics would 
only be identical to the differential equation if At were very small. However, for our purposes 
the continuous and discrete time versions provide an equally good description of an ecosystem 
and we do not wish to distinguish between them. A key point about the equations is that 
the stationary values of the population sizes from equation (|T^) are identical to those from 
equation (|]). 

In the program it is necessary to continuously update the efforts for each species so that 
they remain close to the ESS values. We assume that efforts can change on a time scale of days, 
which is much quicker than the change in population sizes, which occurs on the time scale of 
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the generation time of the organism. If the efforts satisfy Eq. (0), then they also satisfy the 
ESS condition that the gain per unit effort is equal for all prey. If we begin with some choice of 
efforts that is not the ESS and substitute these into the Qij functions on the right of equation 
( |T4|) we obtain a new set of efforts that is closer to the ESS. Repeated iteration of this equation 
therefore causes the efforts to converge on the ESS. In our simulations, after each update of 
the population sizes using equation (]IB|), we updated the efforts by iterating equations (|I^) 



and many times whilst keeping the population sizes fixed, until the maximum relative 
change in effort was smaller than a threshold. In most of our simulations, this threshold was 
10%. Only then did we proceed with the next update of the population sizes. 

In principle, a species i can assign part of its effort to any species j for which Sij > 0, 
however, the ESS condition means that not all such species end up with a non-zero fraction 



of the effort. From ([TSD , the gain per unit effort has a limit as fij tends to 0, and this is the 
maximum achievable value of Qij/fij- If this maximum value is less than the gain per unit 
effort that can be achieved from some other combination of prey excluding j, then the ESS 
solution has fij = 0, i.e. i does not include j in its diet. We believe from numerical studies 
that there is a unique ESS choice of efforts for any fixed set of population sizes. We have 
proved this in some special cases, but have not yet found a general proof. 

The efforts of each species change continuously during the simulation. If = at some 
point in time, it does not necessarily remain so. If the population size of species j increases 
it may pay i to switch some of its effort to j. Also if a third species goes extinct which was a 
well-adapted predator of j, it may be possible for i to feed on j, whereas previously it could 
not do so because it was outcompeted by the third species. In cases where the ESS solution 



for a particular effort is zero, iteration of equation (|T^), from a small starting value, causes it 



to become ever closer to zero, and eventually the computer sets it to zero when it falls below 
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the smallest real number allowed by the operating system. This creates a problem, since if an 
effort is exactly zero, it can never increase again by iteration of equation (Q). Therefore we 
introduced a minimum effort /min in the simulation program, such that whenever the value of 
an effort fij became smaller than /min, but the score Sij was positive, we set fij = fmm- This 
allows efforts that were previously effectively zero to recover again to large values if conditions 
change. In cases where the score Sij is zero, then the corresponding effort is also exactly zero. 
In these simulations we chose /min = 10~^, and found that the results do not depend on the 
precise value of /min, as long as it is small enough. 

Beginning from the same ratio-dependent functional response for one predator and one 
prey (equation Michalski & Arditi (1995) and Arditi & Michalski (1995) generalised the 
equation to a food web in a way that is different from our equation (|15]). These authors 
introduced quantities Xl^''\ the part of species i that is being accessed as a resource by 
species j, and Xj^^\ the part of species j that is acting as a consumer of species i. Similar 
quantities were used by Berryman et al. (1995) with a different form of the functional response. 
The values of these quantities must be determined self-consistently at each moment in time. 
This is qualitatively similar to the way in which we determine the solution for the efforts 
fij at each moment in time. Our equations are simpler because they only require one set of 
auxiliary variables per species rather than two. We are also able to give an interpretation of 
the fij in terms of the ESS, which was not done with the alternative formulation using the 
X parameters. The two formulations nevertheless predict similar effects. Michalski & Arditi 
(1995) show that as the X parameters change, links can appear and disappear from the food 
web, and hence that the structure of the web is not the same at equilibrium as away from 
equilibrium. The same behaviour is seen in our approach, since the fij can change from zero 
to non-zero and vice-versa. 
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5 Competition for external resources in absence of pre- 
dation 



We begin by considering the competition between primary producer species for the external 
resources R in the absence of consumer species at higher levels of the food web. We wish to 
determine how great a degree of diversity can be generated by evolution in this case, and to 
study the way the competition strength affects this diversity. All species feed on the external 
resources ("species 0") only, therefore SiQ > 0. All the scores Sij for interaction between 
species are set to zero in this case. Since each species has only one food source, all the efforts 
are identical to 1. We initialize the model by assigning the values of the feature score matrix 
maf3, and by choosing 10 random features to be the features of the environment. We then 
introduce the first species with 10 randomly chosen features and a population size A'^i(O) = 1. 
After iterating the population equations until all population sizes converge, a new species 
is created by speciation, as described in section 0, and this process is repeated for many 
evolutionary time steps. 

For survival as a primary producer, a species must have a score Sio > b/X, as shown by 
equation above. In addition species can coexist only if their scores are sufficiently close. 
The conditions for coexistence in the case of just two species were given in equation ([10|). 
Species with scores which are too low are out-competed, and become extinct. The strength 
of competition aij between species depends on their degree of similarity qij as defined in 
equation (|l3)- Thus a species with a relatively low score that is phenotypically distant from its 
competitors (shares few features) experiences reduced competition and may survive, whereas 
another species with the same score that is similar to a well-adapted high-score species may be 
out-competed. The rationale behind this is that different species can use resources in different 
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ways. If plants diversify by being of different sizes, by adapting to different temperatures 
and moisture levels, and by adopting different means of dispersal etc., they can make more 
effective use of the fixed amount of sunlight and ground space that is available. When c is 
close to 1 there is strong competition even between distantly relates species. When c is small 
there is much weaker competition between distant species, hence we would expect greater 
diversification of the ecosystem when c is small. 

The results of one simulation run are shown in Figure 1. One can see that the species 
configuration becomes fixed after approximately 13000 time steps. We continued the simula- 
tion to 100000 time steps and did not see any change. Obviously, the species configuration is 
such that no new species can be generated that can survive in the presence of all the existing 
species. However, the new species generated are all similar to existing species (they differ by 
only one feature out of 10 from their parent species). Thus the set of species that arises is 
not necessarily stable against all possible species, just against those which can arise by small 
changes of the existing set. If we use the same score matrix rriais and the same set of fea- 
tures for the external environment, but we start with a different initial species, the stationary 
configuration is different. We also found that the surviving species with the highest score is 
usually not the one with the largest population. This is because of the dependence of com- 
petition on similarity. The species with the highest possible score against the environment is 
usually not part of the stationary configuration. For example, when c = 0, we found that the 
mean score of the species in the stationary configuration is below 6.0, while the best possible 
score is close to 8.0 for our choice of the environment and the score matrix. The similarity 
between the majority of pairs of species was = 0. When the simulation was started with 
the species with the highest score, it died out quickly, but the mean stationary score was 
higher than with a random initial species. This situation was different when c was larger. For 
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c > 0.8, the advantage of being different was smaller, and the stationary population contained 
the species with the highest possible score, together with other species with a high degree of 
similarity to it. 

These results in absence of predation confirm the importance of the parameter c, and 
show that reduction in competition as similarity between species decreases is an important 
factor in promoting the evolution of diversity. They also show that in absence of predation 
the ecosystem evolves towards frozen state that cannot be invaded by new species. 

6 Webs with predation 

In this section we simulate the full Webworld model in the presence of predation. We start the 
simulations again with a single primary producer species, but as the species diversify, webs 
with several trophic layers are created. Even though the network of interactions between 
species can be complex, we find that the iteration of the population dynamics equations 
usually converges rapidly to a fixed point. This property was also shared by the simpler set 
we used in Paper I. This means that we can wait for convergence of all the population sizes for 
each set of species before the next species is added. More complex dynamics is found in other 
food web models such as the one discussed by Blasius et al (1999). Our main interest here 
is to determine which species survive in the in the long-term. Thus the important features 
of the dynamics are whether the newly evolved species can increase in number initially, and 
whether any other species drop below the extinction threshold, N — 1. 

The following figures and data show a selection of simulation results. Figure 2 shows the 
number of species in the system as function of time, measured again in speciation events. The 
two curves differ in the set of random numbers, but not in the parameter values. This means 
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that they differ in their score matrix, in the features of the external environment and of the 
first species. 

In contrast to the case without predation, the web has a continuous overturn of species 
even after a long time. This result is different from Paper I, where the species configuration 
became so well adapted that it could not be invaded by new species. From Fig. ^ one can 
also see that simulations with different random numbers give rise to webs with similar species 
numbers and fluctuation strength. It appears that the differences between runs with the 
same parameter values are smaller than they were with the previous equations (c.f. figure 3 
in Paper I). We also looked at other quantities besides the species numbers, such as those 
shown in Table 1 below, and they are similar for the two simulation runs. 

Next, we studied the influence of the parameter c in Eq. ( ]TlD on the properties of the 
web. As discussed earlier in this paper, a smaller value of c leads to less competition between 
species, and it promotes diversity. Figure |^ shows the number of species as function of time 
for four different values of the parameter c, all other parameters being equal. One can see 
that the species number increases with decreasing c, due to the decrease in competition. We 
have argued before that as c decreases, the efficiency in exploiting a food source depends more 
on the overlap of a species with its competitors, and less on its score. This is demonstrated 
in Figure which shows selected scores as function of time. Initially, the basal species with 
the highest score Sio is chosen and its score is plotted as long as the species exists. When 
the monitored basal species becomes extinct, the basal species with the best score at that 
moment is chosen and monitored, and so on. We have done the same for the predator-prey 
pair with the highest score Sij. Each step in the curves means that the monitored species has 
become extinct. The figure also shows that species overturn is higher on higher trophic levels. 
This is no surprise, since basal species have the largest population sizes, and it is therefore 
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more difficult to drive them to extinction. We find that the scores are higher for larger c, and 
that, in particular, the basal species are replaced less often for larger c. For a value of c = 0.8 
there was an indication that the monitored basal species had become fixed. 

In Paper I we looked extensively at the structure and statistical properties of the food 
webs generated by Webworld, and compared these to real food webs. We now wish to look at 
these properties in the new simulations to see how the changes in the form of the dynamical 
equations influence the web structure. It turns out that, as in Paper I, very reasonable 
agreement with quantities measured in real food webs can be achieved for certain values of 
the parameters of the model. For the purpose of defining food web structure, we consider a 
link between species i and species j to be present if species i consumes at least one individual 
of species j per unit time, i.e. if Qij > 1. As in Paper I, we define the trophic level of a 
species to be the number of links on the shortest path from the external resources to that 
species. In the results tables the Average level is the mean trophic level averaged over all 
species in the web and averaged over time. The Average maximum level is the mean value 
of the maximum trophic level in the web averaged over time. In the analysis of food webs, 
species are often classified according to whether they are basal, intermediate, or top. Basal 
species live exclusively on the external resources (i.e. they have no prey). Top species have no 
predators. Intermediate species have both predators and prey. The following table summarizes 
the results for different values of c, and for the same parameters as in Fig. |^, averaged over 
several thousand time steps: 
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c 
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0.4 


0.2 


No. of species 


27 




55 




79 
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Links per species 


1.68 


1 


70 
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33 
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Av. level 
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69 
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Basal species (%) 


12 




Q 




8 

(J 
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Intermediate species (%) 


86 




90 




90 


90 


Top species (%) 


2 




1 




2 


8 


Mean overlap level 1 


0.71 
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22 


0.06 


Mean overlap level 2 


0.30 





13 





08 


0.04 


Mean overlap level 3 


0.17 





09 





07 


0.04 



Table 1. Results of simulations of the model with R = 10^ and 6 = 5 x 10 ^ 
for four values of the competition parameter c. 

As can be seen from Fig. ^, the two simulations for c = 0.4 and c = 0.2 have not yet reached 
their stationary state. Nevertheless, Table 1 shows several trends present: with decreasing c, 
the fraction of intermediate species, the number of links per species, and the average trophic 
level of a species increases, whilst the fraction of basal species decreases. The mean overlaps 
on levels 1, 2 and 3 are the mean values of the quantity qij (fraction of shared features) for 
all pairs of species on the same level. We observe that the overlap is higher on the lower 
levels (i.e. lower level species are more diverse). The mean overlap on each level decreases as 
c decreases, because the strength of competition increases. The same effect was discussed in 
the previous section for the case with no predation. 

The effect of the size of the resources on the number of species is shown in Figure |^. As 
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would be expected, a larger set of resources can sustain a larger web of species. Table 2 shows 
the mean values of selected properties of the web for these simulations: 
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1.0 X 10^ 1.0 X 10^ 3.5 X 10^ 1.0 X 10*^ 


iNO. oi species 
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1.91 
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Av. level 
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Av. max. level 
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Intermediate species (%) 
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89 


89 


Top species (%) 


2 


2 


6 


1 


Mean overlap level 1 


0.32 


0.34 


0.31 


0.27 


Mean overlap level 2 


0.17 


0.12 


0.11 


0.15 


Mean overlap level 3 


0.19 


0.09 


0.09 


0.12 



Table 2. Results of simulations of the model with c = 0.5 and 6 = 5 x 10 



for four values of the resource R. 

With increasing size of resources, the number of species, the number of levels, and the 
number of links per species increase. A larger fraction of species are intermediate species. An 
exception is the last simulation (with 10® resources) , which has not yet reached its stationary 
state. 

We also studied the dependence of the model properties on the other parameters. As b is 
increased, it is more difficult for a species to become established, and in particular during the 
early stages of a simulation run, the species numbers are smaller. 

When the expression for the aij in Eq. (|Tl]) is modified, the simulation results are qual- 

24 



itatively similar. We tested explicitly the choice where aij is 1 for qij = 1, and a constant 
smaller than 1 for all other qij. 

The figures given in the tables represent averages over several runs. The standard deviation 
is moderate for links per species and the average level and average maximum level (about 
5% of the mean), but larger for the number of species and the overlaps (about 10%). The 
fluctuations in the fraction of basal and intermediate species are 10% and 5% respectively, 
similar to what was found in I. Not surprisingly, given the small numbers involved, the top 
species have a large standard deviation (about 50% of the mean). For the simulation with 
the largest value of R, which has not yet fully reached the stationary state, and where the 
fluctuations in the number of species are rather large, the figures given in the table are very 
rough. In Paper I we made extensive comparisons with the statistics of real food webs, hence 
we will not do that here. Similar food webs are generated using the equations in the present 
paper to those in paper I. 

As mentioned in Section 2, the rule for speciation can be chosen in different ways. Instead 
of choosing each species with the same probability to be the parent of a new species, we also 
did a simulation, where a species was chosen with a probability proportional to its population 
size. The mean number of species is smaller if species undergo speciation in proportion to 
their population size. The reason is that more change is happening on the lower trophic levels, 
making it more difficult for species in the higher trophic levels to become established. 

Finally, we studied the size distribution of extinction events. Figure ^ shows the number 
N{s) of events for which s species went extinct during one time step for one long simulation 
run. There is a sharp maximum at s = 1, which is due to the fact that more than 90% 
of the species created by a mutation cannot survive in the presence of all the other species. 
The curve has an exponential decay for larger s, indicating that large extinction events are 
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unlikely. This is very different from the "self-organised critical" behaviour found by Bak and 
Sneppen (1993) in computer simulations of a much simpler model for large-scale evolution, 
where the size distribution of extinction events follows a power law N{s) ~ s~'^ with r ~ 2. 

7 Conclusion 

In this paper, we have studied a model for evolving food webs. We have established a set of 
coupled ecological equations for the population sizes of the different species in a web which 
satisfies the logical requirements put forward by Arditi & Michalski (1995), and in which the 
distribution of foraging effort for each predator follows an evolutionarily stable strategy. We 
have shown that the model generates food web structures that are comparable to those of real 
webs, and have considered the trends in the web statistics with changing parameter values. 

In the absence of predation, the models gives rise to a stable set of species that cannot 
be invaded by any close variant species. In contrast, in the presence of predation, a web 
is built that has a continuous overturn of species. This result is different from the stable 
species configurations found in Paper I. As the size distribution of extinction events falls off 
exponentially, our results are also different from those of several simpler models for large-scale 
evolution, which usually have a size distribution of extinction events that falls of like a power 
law with an exponent close to 2 (Bak & Sneppen, 1993; Amaral & Meyer, 1999). 

A central point in the theory of self-organised critical systems is that small perturbations 
can sometimes lead to large responses. In the original sandpile model (Bak et al, 1988) 
the addition of one sand grain to the top of a pile can sometimes lead to avalanches of 
falling grains. In Webworld, the equivalent effect is the addition of a new species which can 
occasionally lead to several other species becoming extinct. It is important that extinctions 
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occur as a result of the changes in population sizes caused by adding the new species. There 
are no random extinctions in Webworld: we do not remove a species unless its population 
falls below 1. There is also no random replacement of species; when a new species is added, 
the parent species is not removed. It may happen that the new species out-competes the 
parent species and thus replaces it, but this only happens if the new species is better adapted 
than the old one. In contrast, most other macroevolution models (e.g. Bak & Sneppen, 
1993; Sole et al. 1996; Amaral & Meyer, 1999) either include random extinction or random 
replacement of species. If this is done then sooner or later a very well adapted species with 
a high population will be removed by chance, and this is likely to have a large effect on the 
structure of the ecosystem and maybe lead to further extinctions. In the sandpile analogy 
this is like removing a grain from the bottom of the pile. It would not be surprising if changes 
of this type caused large avalanches. It could be argued that chance extinctions might occur 
due to stochastic fluctuations in population sizes. Our dynamics is deterministic and thus 
does not allow for this possibility, however stochastic fluctuations are unlikely to affect species 
with high population sizes sufficiently to drive them extinct. Therefore we feel that simply 
removing species whose populations fall below the threshold value of 1 is an adequate way of 
dealing with extinctions. 

One of the major questions that one would wish to address with models such as ours 
is whether the large scale extinction events observed in the fossil record could arise as a 
result of the internal dynamics of the ecosystem, or whether external causes are required. 
We have found in preliminary simulations with our model (results not shown) that even 
relatively minor changes to the external environment (species in our model) are capable 
of causing large scale extinction events, which in turn lead to the potential for the growth 
of new species. Thus it seems clear to us that external perturbations can cause extinction 
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avalanches. The more interesting question is therefore whether mass extinctions occur with 
a static external environment - the case considered in this paper. Although the long-term 
behaviour of our model shows a continual overturn of species, no large scale extinctions are 
seen and we therefore deduce that environmental changes are required to produce these. As 
we have stressed above, we believe that what is referred to as "internal dynamics" in some 
models is effectively external, since in these models perfectly well adapted species are removed 
by random changes. So in our view these effects are indistinguishable from the elimination of 
species due to some random external perturbation due to environmental change. We conclude 
that great care has to be taken to distinguish between truly internal dynamics and external 
influences. For this a realistic model of evolutionary dynamics is required. It is likely that 
both internal and external effects exist in the spectrum of extinction events seen in the fossil 
record. 

Another question of interest concerns the robustness of the simulation results to modifica- 
tions of the model. As mentioned at several places in this paper, we found that our qualitative 
results are insensitive to a variety of changes that we made. However, we should also mention 
that our findings depend sensitively on a good implementation of the rule for updating the 
efforts. If the efforts are not given enough time to equilibrate with respect to the population 
sizes, small values of the efforts cannot recover quickly when a prey becomes more abundant, 
and we found that this occasionally led to large extinction avalanches which destroyed almost 
the entire web. Another type of undesirable behaviour was also observed in some simulations 
where we did not allow the efforts to equilibrate properly. If the efforts of new species are 
initialized such that they are not close to their equilibrium value, the species configuration 
becomes frozen after some time, because no new species can become established. 

There are many other questions to ask: experimental ones which relate to the comparison 
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with real systems and theoretical ones which have to do with model structure. We hope to 
investigate such questions in the future. However, we believe that the model introduced in 
this paper is an important step in our understanding of the evolution of food web structure, 
being both simple enough to give an understanding of the basic mechanisms at work and 
realistic enough to allow comparison with data collected by ecologists. 
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Appendix. Evolutionarily Stable Strategies 

Here we consider a predator species i and we determine the ESS choice of efforts. Let the total 
population be N^, and suppose that the majority of the population, — rii, have a foraging 
strategy defined by the efforts fij, whilst a small minority, rij, have a different strategy hij. 
Following the usual argument of evolutionary game theory, we require to calculate the payoff 
to the minority and majority strategies, and hence to determine conditions under which the 
minority can invade. In this case the payoff is the total rate of gain of resources from all prey. 
The payoff for the strategy fij in absence of the minority strategy is 
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where, for convenience, we use Kij to denote all the terms in the denominator of the gij 
function in equation (|15D that do not depend on the efforts of species i: 

K,^ = bNj + ^ ak^SkJfkjNk. ( 18) 

The payoff for the minority species in the presence of the majority is: 

j Sij fiji^Ni Tij) + SijhijTii + Kij 

In the above equation, since the two strategies are played by different individuals which are 
members of the same species, the a factor for competition between individuals with different 
strategies is 1, as it is for individuals using the same strategy. We require the payoff to the 
invading strategy when it is rare (i.e. when -C Ni), which is just obtained by setting 
equal to zero in (|19|). In a similar way the payoff to the majority strategy in the presence of 
the minority can be written down, but this reduces to ( |T7| ) when <C A^j. This gives 

Gmm - G = E ^^(^i-/^^) = Y.ih., _ (20) 

j >->ljJlj^^l "T J^lj j Jlj 

where Qij/fij in the equation above is the gain per unit effort from prey j for the majority 
strategy. The minority strategy can invade if Gmin — G > 0. 

Now suppose that the invading strategy differs from the majority strategy for two prey 
species k and I, so that hik = fik + A/, hu = fa — A/, and hij = fij for all the other prey j. 
In this case 

Grmn -G = Af(f-f-]. (21) 
\ Jik Jil ) 

If the gain per unit effort from prey k is higher than that from prey / then any strategy with 
positive A/ can invade, whilst if the reverse is true then any strategy with negative A/ can 
invade. However, if the gain per unit effort is equal for the two prey then variant strategies 
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are neutral. The same can be said for any pair of prey species k and /. It therefore follows 
that the ESS is the strategy with the gain per unit effort being equal for all prey. If neutral 
variant strategies accumulate to a non-negligible fraction, then selection will again operate to 
drive the population back to the ESS. 

It is interesting to note that the ESS does not correspond to the solution predicted by 
optimal foraging theory (OFT) (Stephens & Krebs, 1986). The OFT solution would be to 
maximise G in equation (|T^ with the constraint that the efforts sum to 1. This can be 
calculated, and the result is (by definition) greater than the total gain to a predator when all 
adopt the ESS. However, a single ESS predator in a population of OFT predators actually 
has a higher total gain than the OFT population. Thus the ESS can invade the OFT solution, 
but the reverse is not true. Hence we argue that the ESS is the appropriate choice of efforts 
for our model. 

In most of the models considered by Stephens & Krebs (1986) the payoff to the predator 
is not affected by what other predators do, therefore the straightforward OFT solution of 
optimising the total rate of energy intake is appropriate (see their comments on p 211 regarding 
game theory). However, competition between predators of the same species and between 
different species is an essential part of the way our population dynamics equations are set up, 
and we also believe it is an important factor in real ecosystems. Therefore it is important 
to treat the foraging problem from a game theory point of view. The need for game theory 
has also been defended recently by Reeve & Dugatkin (1998). Various game theory models 
dealing with aspects of foraging behaviour have been proposed (Matsuda et al., 1996; Shaw 
et al., 1995; Leonardsson & Johansson, 1997; Visser, 1991; Giraldeau & Livoreil, 1998; Sih, 
1998). 
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Figure 1: Number of species (solid line) and maximum score (dotted), as function of time for 
R — 1000, b — 0.1, c — 0.5. Time is measured in units of number of speciation events. 
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Figure 2: Number of species as function of time for R = 350000, b = 0.005, c — 0.5. Time is 
measured in units of mutations. The two simulations were carried out with different sets of 
random numbers. 
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Figure 3: Number of species as function of time for R = 100000, b = 0.005, and 
c = 0.2, 0.4, 0.6, 0.8 (from top to bottom curve). Time is measured in units of mutations. 
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Figure 4: Score of the monitored basal species (solid line) and predator-prey pair (dotted 
line) as function of time. The parameters are R — 100000, b — 0.005, and c = 0.5. 
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Figure 5: Number of species as function of time for c = 0.5, b — 0.005, and R — 10^, 10^, 3.5 x 
10^, 10^ (from bottom to top curve). Time is measured in units of mutations. 
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Figure 6: Size distribution of extinction events for a simulation with R — 100000, c = 0.5, 
b = 0.005. Species were chosen with a probabihty proportional to their population size to be 
the parent of a new species. 



40 



